home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / ARGONET / PD / MATHS / RLAB / RLAB125.ZIP / !RLaB / toolbox / autospec < prev    next >
Text File  |  1995-02-25  |  2KB  |  98 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Compute Auto-Spectra of the input sequence
  4.  
  5. // Syntax:    autospec ( X, dt, M, WIN, raw );
  6.  
  7. // Description:
  8.  
  9. //      Compute the auto-spectra (psd) of the input vector X. autospec
  10. //      returns a two-column matrix. The 1st column is the frequency
  11. //      scale and the second column is the values of the one-sided
  12. //      spectra.
  13.  
  14. //    X:    Input sequence
  15. //    dt:    Sample interval 
  16. //              (default = 1)
  17. //    M:    Length of sequence to use for computations. The
  18. //        sequence X is "cut-up" into M-point sequences
  19. //              (default = length (X))
  20. //      WIN:    Window type: "hann", "rect", "hamm"
  21. //              (default = "rect")
  22. //      raw:    If raw == 1, then output the "raw" auto-spectra.
  23. //              Otherwise output the auto-spectra from 0-Nyq/2
  24. //              with a time-scale.
  25. //              (default = 0)
  26.  
  27. // Dependencies
  28.  
  29.    rfile window
  30.    rfile detrend
  31.  
  32. //-------------------------------------------------------------------//
  33.  
  34. autospec = function ( X , dt, M, WIN, raw )
  35. {
  36.   local ( X , dt, M, WIN, raw )
  37.  
  38.   X = X[:];        // Make sure it is a row
  39.   if (!isreal (X)) { error ("autospec: X must be real"); }
  40.   n = length (X);
  41.  
  42.   //
  43.   // Set the defaults.
  44.   //
  45.  
  46.   if (!exist (dt)) { dt = 1; }
  47.   if (!exist (M)) { M = n; }
  48.   if (!exist (WIN)) { WIN = "rect"; }
  49.   if (!exist (raw)) { raw = 0; }
  50.  
  51.   // Some error checking
  52.  
  53.   if (M > n) { error ("spectrum: M > n not allowed"); }
  54.   if (M < 0) { error ("spectrum: M < 0 not allowed"); }
  55.  
  56.   // Create the frequncy scale (in Hertz) up to the Nyquist
  57.  
  58.   fscale = (((0:M-1)/M)/dt)';
  59.  
  60.   k = fix (n/M);    // The number of windows
  61.   index = 1:M;
  62.  
  63.   // Create the window
  64.  
  65.   w = window (M, WIN);          // Window specification
  66.   U = k*sqrt(sum(w.^2)/w.n);    // Normalizing scale factor
  67.  
  68.   Pxx = zeros (M, 1);
  69.  
  70.   for (i in 1:k)
  71.   {
  72.     xw = w .* detrend (X[index]);
  73.     index = index + M;
  74.     Xx = (2/(M*dt)) * abs (dt * fft (xw)) .^ 2;
  75.     Pxx = Pxx + Xx;
  76.   }
  77.  
  78.   if (raw)
  79.   {
  80.  
  81.     //
  82.     // Return the "raw" spectrum.
  83.     //
  84.  
  85.     return Pxx/U;
  86.   else
  87.  
  88.     //
  89.     // Select 1st half and eliminate DC value.
  90.     // Return a matrix with the frequency scale in the 1st column
  91.     // and the PSD in the 2nd.
  92.     //
  93.  
  94.     Pxx = Pxx[2:(M/2)+1]/U;
  95.     return [fscale[2:(M/2)+1], Pxx];
  96.   }
  97. };
  98.